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Performance  of  the  adjoint  and  adjoint-free  4-dimensional  variational  (4dVar)  data  assimilation  techniques 
is  compared  in  application  to  the  hydrographic  surveys  and  velocity  observations  collected  in  the  Adriatic 
Sea  in  2006.  Assimilating  the  data  into  the  Navy  Coastal  Ocean  Model  (NCOM)  has  shown  that  both  methods 
deliver  similar  reduction  of  the  cost  function  and  demonstrate  comparable  forecast  skill  at  approximately  the 
same  computational  expense.  The  obtained  optimal  states  were,  however,  significantly  different  in  terms  of 
distance  from  the  background  state:  application  of  the  adjoint  method  resulted  in  a  30-40%  larger  departure, 
mostly  due  to  the  excessive  level  of  ageostrophic  motions  in  the  southern  basin  of  the  Sea  that  was  not 
covered  by  observations. 
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1.  Introduction 

In  recent  years  the  ensemble  approach  to  data  assimilation  has 
developed  rapidly  due  the  growth  of  computer  power  made  available 
through  massive  parallelization.  An  attractive  feature  of  the  ensemble 
technique  is  its  ability  to  probe  the  structure  of  a  dynamical  system 
and  assess  the  sensitivity  of  its  outputs  (e.g.,  measured  quantities)  to 
variations  in  poorly  known  inputs  (e.g.,  initial  conditions)  without  us¬ 
ing  the  adjoint  code.  In  particular,  feasibility  of  the  ensemble  method 
to  the  estimation  of  sensitivities  was  demonstrated  in  meteorological 
(Ancell  and  Hakim,  2007;  Torn  and  Hakim,  2008)  and  oceanographic 
(Yaremchuk  and  Martin,  2014)  applications. 

Another  important  advantage  of  the  ensemble  approach  is  the 
opportunity  it  offers  to  treat  the  numerical  model  as  a  black  box, 
thus  avoiding  the  burden  of  the  development  and  maintenance  of  the 
tangent  linear  and  adjoint  codes  required  by  the  variational  meth¬ 
ods  (e.g.,  Le  Dimet  and  Tlagrand,  1986).  Employing  this  property, 
Anderson  et  al.  (2009)  and  Hoteit  et  al.  (2013)  developed  the  data  as¬ 
similation  research  testbed  (DART)  system  on  the  basis  of  the  widely 
used  ensemble  Kalman  filter  (EnKF). 

Recently,  significant  progress  also  has  been  made  in  extending 
the  EnKF  technique  into  the  particle  filtering  framework  (e.g.,  Hoteit 
et  al.,  2012)  and  in  coupling  EnKF  techniques  with  3d  and  4d  vari¬ 
ational  methods  (e.g.,  Zupanski,  2005;  Liu  et  al.,  2008;  Zhang  et  al., 
2009).  Of  particular  interest  in  the  present  context  has  been  the 


*  Corresponding  author.  Tel.:  +1  2286885259. 

E-mail  address:  max.yaremchuk@nrlssc.navy.mil  (M.  Yaremchuk). 

http://dx.doi.Org/10.1016/j.ocemod.2015.10.010 
1463-5003/Published  by  Elsevier  Ltd. 


development  of  the  Maximum  Likelihood  Ensemble  Filter  (Zupanski, 
2005)  based  on  the  explicit  computation  of  the  square  root  of  the 
Hessian  matrix  in  the  subspace  spanned  by  the  ensemble  members. 

Merging  ensemble  approaches  with  variational  techniques  has 
developed  along  two  lines:  1)  improvement  of  the  background  error 
covariances  (BECs)  through  introduction  of  the  ensemble-based  esti¬ 
mates  and/or  their  hybrid  generalizations  (Clayton  et  al.,  2013;  Kuhl 
et  al.,  2013)  and  2)  searching  for  the  optimal  solution  within  the  sub¬ 
spaces  spanned  by  the  leading  error  modes  of  the  BECs,  a  technique 
pursued  by  many  authors  in  the  last  decade  (e.g.,  Liu  et  al.,  2008; 
Zhang  et  al.,  2009;  Zhang  and  Zhang,  2012;  Trevisan  et  al.,  2010). 
This  assumption  implicitly  assumes  that  the  BEC  structure  is  well  de¬ 
scribed  by  these  (possibly  localized)  modes.  More  recently,  perfor¬ 
mance  the  adjoint-free  family  of  methods  (4dEnVar)  based  on  the  for¬ 
mulation  by  Liu  et  al.  (2008)  has  been  compared  with  the  4dVar  tech¬ 
nique  in  the  framework  of  idealized  experiments  with  the  Lorenz-05 
model  (Fairbairn  et  al.,  2014).  The  results  show  a  significantly  better 
performance  of  the  4dEnvar  for  moderate-length  assimilation  win¬ 
dows  with  low-density  observations.  Desroziers  et  al.  (2014)  demon¬ 
strated  a  close  relationship  between  the  4dEnVar  and  4dVar  state 
space  formulations  and  compared  various  implementations  of  4dEn- 
Var  with  4dVar  in  an  idealized  setting. 

The  above  cited  developments  mostly  deal  with  meteorological 
applications,  where  the  ensembles  are  supported  by  significantly 
denser  data  than  are  available  in  the  ocean.  High  data  density  al¬ 
lows  one  to  obtain  reasonably  good  estimates  of  the  BECs  from  the 
ensemble  using  truncated  representation  of  the  localization  matrices 
and  to  efficiently  compute  the  cost  function  gradient  on  the  model 
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grid  directly  from  ensemble  perturbations  (Liu  et  al.,  2009;  Tian  and 
Xie,  2012).  A  significant  advantage  of  such  an  approach  is  the  ab¬ 
sence  of  the  necessity  to  develop  and  maintain  tangent  linear  and 
adjoint  codes  and  its  flexibility  in  adaptation  to  various  dynamical 
constraints. 

In  the  ocean,  the  ensemble-based  BEC  estimates  tend  to  be  less 
accurate,  and  one  has  to  rely  on  ad  hoc  BEC  representations  (e.g., 
Mirouze  and  Weaver,  2010;  Yaremchuk  and  Sentchev,  2012).  De¬ 
velopment  of  an  efficient  adjoint-free  assimilation  method  also  be¬ 
comes  more  problematic  as  one  has  to  select  a  few  reliable  pertur¬ 
bations  with  more  care.  Early  attempts  to  develop  practical  a4dVar 
algorithms  in  oceanography  were  limited  to  predetermined  low¬ 
dimensional  subspaces  spanned  either  by  the  reduced-order  approx¬ 
imations  of  the  model  Green’s  functions  (Stammer  and  Wunsch, 
1996;  Menemenlis  and  Wunsch,  1997),  or  by  the  dominant  princi¬ 
pal  component  vectors  (EOFs)  associated  with  the  model  statistics 
(e.g.,  Robert  et  al.,  2005;  Qui  et  al.,  2007;  Hoteit,  2008).  In  fact,  the 
4dEnVar  technique  pursues  a  similar,  but  more  general  approach,  pa¬ 
rameterizing  the  search  subspace  by  Schur  products  of  the  ensemble 
members  with  the  eigenvectors  of  the  reduced-order  representation 
of  the  localization  matrix. 

In  the  present  paper,  the  a4dVar  approach  of  Yaremchuk  et  al. 
(2009)  is  tested  against  the  observation  space  4dVar.  Both  algo¬ 
rithms  are  dynamically  constrained  by  the  Navy  Coastal  Ocean  Model 
(NCOM)  and  applied  to  the  real  observations  collected  in  the  Adriatic 
Sea  in  August,  2006.  A  specific  feature  of  the  a4dVar  method  tested 
here  is  that  it  employs  an  iterative  search  over  a  sequence  of  low¬ 
dimensional  subspaces  to  find  the  cost  function  minimum. 

The  paper  is  organized  as  follows.  In  the  next  section  we  briefly 
describe  the  4dVar  methodology,  outline  the  a4dVar  method  to  be 
tested,  and  provide  a  detailed  description  of  the  experimental  set¬ 
ting.  In  Section  3  performance  of  the  4dVar  and  a4dVar  methods  is 
compared  in  terms  of  the  convergence  rate,  forecast  skill,  computa¬ 
tional  expense,  and  particular  properties  of  the  optimized  solutions. 
Conclusions  and  discussion  of  the  results  are  presented  in  Section  4. 


2.  Methodology 

2.1.  4dVar  assimilation 


For  the  sake  of  clarity,  consider  the  4dVar  method  as  the  following 
linear  discrete  least-squares  problem  constrained  by  model  dynamics 
in  a  small  vicinity  of  the  model’s  background  trajectory  xj|: 


rET 'x°  +  £  (H„xn  -  d")^-1  (Hnxn  -  d") 


min.  (1) 


where  n  enumerates  observation  times,  B  is  the  background  error  co- 
variance  matrix  of  x°  which  describes  the  (Gaussian)  error  statistics 
of  the  model  state  at  n  =  0,  H„  is  the  model-data  projection  opera¬ 
tor,  d'1  is  the  misfit  y(J  -  H„xjJ  between  observations  y°n  and  the  cor¬ 
responding  background  model  values,  R„  is  the  observation  error  co- 
variance,  and  T  denotes  transposition.  Further  below  we  denote  the 
dimension  of  the  discretized  model  state  vector  x  by  M  and  the  total 
number  of  observations  by  K. 

The  optimal  correction  vector  xn  is  governed  by  the  recursive  re¬ 
lationship  xn  =  Mnxn_1,  where  M„  is  the  dynamical  operator  of  the 
model  linearized  in  the  vicinity  of  the  background  trajectory  xjj  at  the 
time  interval  (t„_1 ,  tn),  so  that 

xn  =  M,1M„_1...M2M1x0.  (2) 

Introducing  new  notation  c  =  x°  for  the  control  vector,  M"  = 
M„  . . .  M2Mt  for  the  aggregated  n-step  propagator,  H„  =  R„  1/2H„, 
d"  =  R,j1/2dn,  omitting  over-bars,  and  taking  (2)  into  account,  the 
minimization  problem  (1)  can  be  rewritten  in  terms  of  the  optimal 


correction  c  to  the  initial  state: 


J=l 


cTB  ’c  +  (H„M"c  -  dn)T(H„Mnc  -  d") 


min.  (3) 


A  4dVar  data  assimilation  method  finds  the  minimum  of ]  by  solv¬ 
ing  the  normal  equation: 


VJ  =  B  :c  +  Y]  MnTHj (H„M"c  -  d")  =  0. 


(4) 


To  simplify  further  treatment,  introduce  the  following  notation  for 
the  Hessian  matrix  H  and  the  right-hand  side  b, 

H  =  B_1  +^M',thJh„M";  b  =  J]MnTH^dn,  (5) 

n  n 


which  define  the  solution  of  the  normal  equation  He  -  b  =  0. 


2.3.1.  Adjoint  techniques 

There  are  two  major  approaches  to  4dVar  assimilation.  The  first 
one  (the  so-called  “state  space  approach”)  typically  solves  (4)  iter¬ 
atively  through  a  conjugate  gradient  descent  or  related  algorithm, 
which  requires  on  every  iteration  the  formation  of  a  matrix-vector 
product  using  the  Hessian,  H,  or  an  equivalent  process,  which  in  ei¬ 
ther  case  involves  application  of  the  non-linear  model  operator,  MT 
(the  “adjoint  model")  and  (in  many  cases)  the  linearized  model  M. 
The  method  is  widely  used  in  a  number  of  community  OGCMs  (MIT, 
ROMS),  and  in  operational  meteorology  (ECMWF). 

Such  iterative  solution  approaches  generate  a  sequence  of  residu¬ 
als  r;  =  He,  -  b.  which  is  the  cost  function  gradient  (4)  evaluated  at 
the  current  solution  iterate,  c,-.  However,  knowledge  of  the  steepest 
descent  direction  for  the  cost  function  requires  access  to  the  adjoint 
model,  Mt;  note  that  it  enters  the  expressions  (5)  for  both  the  Hes¬ 
sian  and  the  right-hand  side  of  the  normal  system. 

Numerically,  the  procedure  of  calculating  r,-  involves  two  major 
steps: 

(1)  Sequential  calculation  of  xj1  =  Mnc>  (forward  run  of  the  tan¬ 
gent  linear  model)  supplemented  by  additional  calculation  of 
the  model-data  misfits 

qj1  =  (Hnx"  —  d").  (6) 

(2)  Summation  of  the  products  M"Tq"  conveniently  performed  in 
reverse-time  order,  because  MnT  =  (M„  . . .  Mi)t  =  M^. . .  M)[. 
This  corresponds  to  backward-in-time  integration  of  the  ad¬ 
joint  model  forced  by  qj1. 

The  second  approach  to  4dVar  (“observation  space  method")  ap¬ 
pears  to  be  more  rigorous  in  that  it  can  include  more  easily  ex¬ 
plicit  treatment  of  the  model  errors  e"  =  xn  -  M^x"-1.  In  this  for¬ 
mulation,  the  cost  function  (1)  is  augmented  with  an  additional 
term  enTBn1en  involving  model  errors.  In  other  words,  the  back¬ 
ground  error  is  explicitly  separated  into  the  components  associated 
with  the  uncertainty  B0  of  the  initial  state,  and  the  uncertainty  B„ 
of  the  model  equations/forcing.  The  normal  equation  in  this  case 
is  more  complicated  than  (4),  and  can  be  solved  numerically  us¬ 
ing  the  representer  method  (e.g.,  Bennett,  2002;  Rosmond  and  Xu, 
2006).  The  latter  is  closely  related  to  the  optimal  interpolation  tech¬ 
nique  as  it  seeks  a  solution  of  the  normal  equation  in  the  form  of  a 
linear  mapping  of  the  model-data  misfits,  dn  on  the  control  space. 
Mathematically,  the  approach  employs  the  Sherman-Morrison- 
Woodbury  formula  to  transform  the  Hessian  inverse  from  the  state 
space  to  the  observational  space,  which  has  a  smaller  dimension 
in  oceanographic  applications.  On  the  other  hand,  minimization  of 
the  (non-linear)  cost  function  in  the  observation  space  contains  two 
embedded  loops,  and  involves  multiple  convolutions  with  the  non- 
sparse  matrix  B,  making  the  method  sometimes  more  computation¬ 
ally  expensive.  It  has  been  implemented  in  ROMS  (Moore  et  al.,  2011 ) 


M.  Yaremchuket aL/ Ocean  Modelling97 (2016)  129-140 


131 


as  an  optional  feature,  and  in  the  Naval  Research  Laboratory  for  both 
atmospheric  (Xu  and  Rosmond,  2004;  Xu  et  al.,  2005)  and  oceanic 
(Ngodock  and  Carrier,  2014)  data  assimilation  systems. 

2.12.  Adjoint-free  methods 

With  the  ongoing  trend  toward  parallelization  in  computer  tech¬ 
nologies,  directly  perturbing  a  large  number  of  control  variables  be¬ 
comes  computationally  feasible,  allowing  the  replacement  of  tangent 
linear  and  adjoint  codes  by  the  finite  differences  of  perturbed  numer¬ 
ical  models  (’’ensemble  members”)  that  may  be  run  in  parallel.  For 
example,  the  fast-developing  4dEnVar  method  (e.g.,  Liu  et  al.,  2009; 
Buehner  et  al,  2010;  Desroziers  et  al.,  2014)  restricts  the  search  for 
an  optimal  increment  to  a  pre-determined  subspace  spanned  by  the 
eigenvectors  of  a  localized  ensemble  covariance  matrix. 

When  skillful  ensembles  are  unavailable,  a  similar  search  could 
be  performed  using  the  leading  eigenvectors  of  an  ad  hoc  model  of 
B,  provided  that  search  directions  (SDs)  are  kept  orthogonal  with 
respect  to  the  inner  product  associated  with  the  Hessian  matrix 
(Appendix  B).  The  adjoint-free  (a4dVar)  method  considered  here  fol¬ 
lows  this  approach  through  iterative  minimization  of  the  cost  func¬ 
tion  in  a  sequence  of  subspaces  spanned  by  the  SDs 

s"  :=  [B"1 +  H^H„]_1q"  =  B[l  +  H^HnB]_1qn;  n  =  0,..,N,  (7) 

where  I  is  the  identity  operator  in  state  space. 

Compared  to  the  search  directions  defined  by  the  leading  eigen¬ 
vectors  of  B,  the  a4dVar  subspaces  spanned  by  {sn}  contain  extra  in¬ 
formation  on  the  structure  of  the  Hessian  (through  the  observation 
operators  H„)  and  on  the  magnitudes  of  model-data  misfits  q"  at  the 
current  iteration.  The  overall  a4dVar  strategy  is  to  replace  the  single 
SD  obtained  by  projecting  model-data  misfits  q"  on  the  initial  condi¬ 
tion  with  the  adjoint  model  by  parallel  searches  in  multiple  directions 
(7).  Although  these  directions  are  unlikely  to  include  the  direction  of 
local  steepest  descent,  they  implicitly  accumulate  information  on  the 
Hessian  structure  through  the  H-orthogonalization  process  and  may 
therefore  be  competitive  in  efficiency  with  the  adjoint  method  in  a 
typical  oceanographic  application. 

Feasibility  of  the  presented  a4dVar  approach  is  also  motivated  by 
the  following  considerations: 

(1)  the  search  subspaces  can  be  easily  H-orthogonalized  to  each 
other  to  avoid  redundant  “nearly  parallel”  searches  in  the  di¬ 
rections  that  have  been  largely  explored,  and,  therefore,  bring 
only  minor  reduction  of  the  cost  function; 

(2)  search  subspaces  spanned  by  {sn}  are  unlikely  to  be  strictly 
H-orthogonal  to  the  cost  function  gradient  and  may  provide 
a  larger  overall  projection  on  the  direction  towards  the  cost 
function  minimum  than  the  direction  of  steepest  descent; 

(3)  adjoint  and  tangent  linear  codes  for  the  state-of-the  art  GCMs 
are  never  exact  and  require  several  times  more  computational 
resources  compared  to  direct  model  runs,  providing  an  incen¬ 
tive  to  avoid  full-adjoint  model  evaluation  and  to  take  advan¬ 
tage  of  the  parallelism  afforded  by  the  multiple  direct  model 
runs  used  by  a4dVar. 

Note  that  the  number  of  SDs  is  not  restricted  to  the  number  of  ob¬ 
servation  times.  In  principle,  a  simultaneous  search  can  be  performed 
in  all  the  directions  corresponding  to  every  single  observation  within 
the  assimilation  window.  However,  care  should  be  taken  to  exclude 
"nearly  parallel"  SDs  from  consideration.  This  can  be  accomplished 
by  extracting  the  leading  eigenvectors  from  the  given  set  of  search 
directions.  In  the  reported  experiments,  we  assigned  a  SD  to  all  ob¬ 
servations  collected  at  an  observation  time. 

Plausible  forms  of  SD  generation  are  not  restricted  to  (7).  For  in¬ 
stance,  search  subspaces  could  be  built  using  the  leading  modes  of 
the  model  trajectory  at  subsequent  iterations  provided  the  SDs  are 
H-orthogonal  and  have  sizable  projections  on  the  local  gradient.  This 


strategy  has  been  widely  used  in  the  early  versions  of  ensemble- 
4dVar  techniques  (Qui  et  al.,  2007;  Hoteit,  2008)  and  proved  to  be 
rather  effective  in  a  number  of  applications  including  the  tested 
a4dVar  method  (Yaremchuk  et  al„  2009;  Panteleev  et  al„  2015). 

As  an  alternative  option,  one  can  simply  use  eigenvectors  of  B  in 
the  descending  order  of  their  eigenvalues  to  build  the  search  sub¬ 
spaces.  In  the  latter  case,  the  convergence  rate  for  the  respective 
a4dVar  algorithm  can  be  estimated  (more  details  on  the  subject  can 
be  found  in  Appendix  A).  However,  this  approach  demonstrated  a  sig¬ 
nificantly  slower  convergence  in  the  realistic  application  considered 
below.  We  partly  attribute  this  behavior  to  the  above  mentioned  ab¬ 
sence  of  information  on  the  (unevenly  distributed)  observational  lo¬ 
cations  in  the  subspaces  spanned  by  the  eigenvectors  of  B. 

The  primary  purpose  of  the  present  study  is  to  test  the  perfor¬ 
mance  of  the  a4dVar  algorithm  in  application  to  a  realistic  oceano¬ 
graphic  data  set  constrained  by  a  state-of-the-art  numerical  model 
and  compare  its  performance  with  a  standard  4dVar  approach.  The 
method  is  based  on  (7)  and  outlined  as  follows: 

0.  Specify  the  dimension  ms  of  the  search  subspaces,  their  num¬ 
ber  k  to  be  kept  in  memory  for  H-orthogonalization,  the  maxi¬ 
mum  number  of  iterations  I,  the  perturbation  magnitude  e  and 
the  background  model  trajectory  xjj.  Set  the  iteration  number  i 
to  zero,  c0  =  0.  and  compute  d". 

1.  Compute  xl Y,-  =  H1/2c,  and  the  search  directions  s?  Eq.  (7). 

2.  Extract  the  ms  leading  EOFs  pf ,  m  =  1 . ms  of  the  search  di¬ 

rections  to  form  the  basis  in  the  search  subspace. 

3.  Perturb  the  initial  conditions  c,  c,  +  epj"  and  run  (in  paral¬ 
lel)  the  ensemble  of  ms  perturbed  models,  computing  the  re¬ 
spective  perturbed  values  of  Jt"  and  Yj". 

3.  H-orthogonalize  the  search  basis  {pi"}  with  respect  to  at  most 
k  basis  vectors  obtained  on  the  previous  iterations  and  com¬ 
pute  optimal  corrections  8ct  (see  Appendix  B). 

4.  Set  ci+1  =  Cj  +  8cj. 

5.  If  i  =  I  exit.  Otherwise  set  i  •«-  i  +  1 ,  then  go  to  1. 

The  above  listed  algorithm  was  configured  for  the  purpose  of  com¬ 
parison  with  4dVar  both  in  terms  of  accuracy  and  computational  ex¬ 
pense.  For  this  reason  the  value  of  I  was  set  by  the  total  CPU  time  re¬ 
quired  for  convergence  of  the  4dVar  algorithm.  Parameters  ms,  k,  and 
e  were  selected  by  experimentation.  The  H-orthogonalization  proce¬ 
dure  and  computation  of  8ct  are  based  on  the  technique  proposed  by 
Zupanski  (2005),  with  more  details  given  in  Appendix  B. 

Comparison  of  4dVar  and  a4dVar  is  done  by  constraining  NCOM 
with  observations  collected  in  the  Adriatic  Sea  in  August,  2006.  De¬ 
scription  of  the  model,  data,  and  parameters  of  the  assimilation  algo¬ 
rithms  are  given  below. 

2.2.  Model  and  data 
2.2.3.  The  model 

The  NCOM  is  a  free-surface  primitive-equation  hydrostatic  ocean 
model  with  a  coordinates  in  the  upper  layers  and,  optionally,  fixed 
depths  below  a  user-specified  distance  from  the  surface.  Algorithms 
that  comprise  a  NCOM  computational  kernel  are  described  in  de¬ 
tail  by  Martin  (2000)  and  with  some  improvements  by  Morey  et  al. 
(2003)  and  Barron  et  al.  (2006).  The  vertical  mixing  model  utilized  is 
the  Mellor-Yamada  level  2  closure  scheme  (Mellor  and  Yamada,  1974) 
and  the  equation  of  state  of  Mellor  (1991)  is  used.  Biharmonic  hor¬ 
izontal  diffusion  is  prescribed  implicitly  via  third-order  upwind  ad- 
vection  (Holland  et  al.,  1998). 

The  model  was  configured  at  3  km  resolution  on  an  85  x  294  hor¬ 
izontal  grid  (Fig.  1)  with  32  levels  in  the  vertical.  The  top  22  a  levels 
follow  the  bathymetry,  stretching  from  the  surface  to  a  fixed  depth  of 
291  m,  and  10  fixed-depth  levels  are  used  below  291  m. 

Initial  and  open  boundary  conditions  for  the  sea  surface  height  £, 
temperature  T,  salinity  S,  and  horizontal  velocities  u.  v  were  provided 
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Fig.  1.  Model  domain  with  CTD  stations  (circles)  and  moorings  (triangles)  of  the  DART 
experiment.  Gray  contours  (m)  show  bottom  topography. 


from  the  global  NCOM  (Barron  et  al„  2004)  solution  for  the  region. 
Tidal  forcing  is  not  used.  The  model  was  forced  by  the  river  runoff  and 
atmospheric  fields  derived  from  the  regional  ALADIN  atmospheric 
model  with  8  km  horizontal  resolution  (Ivatek-Sahdan  and  Tudor, 
2004). 

In  the  described  assimilation  experiments,  initial  conditions  were 
used  as  control  variables,  i.e.,  the  vector  c  comprised  all  the  grid  point 
values  of  f ,  T,  S,  u,  v  at  n  =  0.  With  the  given  3-dimensional  grid  and 
bathymetry,  the  inverse  problem  has  M=l,493,570  unknowns. 

The  first  guess  (background)  values  of  c  were  taken  from  the 
NCOM  simulation  described  by  Martin  et  al.  (2009)  and  then  adjusted 
to  suppress  temperature  and  salinity  biases  during  the  assimilation 
time  interval  (0.00  UTC  on  08/14  to  0.00  UTC  on  08/29/2006).  After 
the  adjustment,  the  horizontal-and-time  average  misfits  between  the 
background  solutions  and  TS  observations  did  not  exceed  0.02  “C  and 
0.005  psu,  respectively. 

Vertical  profiles  of  the  NCOM  standard  deviations  (Fig.  2)  demon¬ 
strate  quantitative  similarity  between  the  variability  of  the  first  guess 
model  solution  and  the  mean  variability  at  the  observation  points. 
However,  the  root-mean  square  (rms)  model-data  misfits  were  found 
to  be  of  the  same  order  of  magnitude  as  the  observed  variability,  indi¬ 
cating  that  the  model  has  limited  simulation  skill  of  the  smaller-scale 
features  without  further  adjustment  of  the  poorly  known  parameters 
(e.g.,  initial  conditions)  provided  by  data  assimilation. 

It  is  also  noteworthy,  that  the  obtained  background  solution  pro¬ 
vides  a  rather  challenging  environment  for  variational  assimilation, 
as  it  contains  a  large  number  of  transient  mesoscale  features  induced 
by  in  wind  forcing  events  and  instabilities  of  the  coastal  boundary 
currents  (see,  for  example,  Figs.  7  and  8  in  Burrage  et  al.,  2009).  In 
particular,  in  the  upper  150  m  (sampled  by  90%  of  observations)  the 


Fig.  2.  Vertical  profiles  of  the  horizontal  and  time  averaged  rms  variations  of  tempera¬ 
ture  and  salinity  (left)  and  velocity  components  (middle)  derived  from  the  background 
model  run  and  observations.  Right  panel  shows  vertical  profiles  of  the  rms  model-data 
differences  normalized  by  the  respective  rms  variabilities  of  observations. 


ratio  of  the  time-averaged  magnitudes  of  non-linear  to  linear  terms 
in  the  momentum  equation  was  close  to  0.3,  indicating  a  consider¬ 
able  degree  of  dependence  of  the  tangent  linear  and  adjoint  models 
on  the  background  solution. 

2.2.2.  DART06  experiment 

The  data  used  in  this  research  were  acquired  in  the  course  of  a 
collaborative  field  experiment  in  the  central  Adriatic,  the  Dynamics 
of  the  Adriatic  in  Real  Time  (DART)  (Martin  et  al.,  2009;  Burrage  et  al., 
2009).  In  the  present  study,  ADCP  and  CTD  observations  from  August 
14  to  August  29,  2006  are  used  (Fig.  1). 

Current  velocities  u.  v  were  measured  by  19  moored  ADCPs  at  lo¬ 
cations  shown  by  triangles  in  Fig.  1.  Due  to  the  inaccuracy  of  ADCP 
measurements  in  proximity  to  the  surface  and  bottom,  the  data 
spanned  the  depth  range  from  15  to  150  m.  All  the  velocity  data  were 
detided  and  averaged  over  29, 12  h  intervals  centered  at  the  assimila¬ 
tion  times  t„  of  0  and  12  UTC.  The  average  duration  of  an  ADCP  time 
series  employed  in  the  data  assimilation  experiments  was  12.7  days. 

Temperature  T  and  salinity  S  were  measured  at  219  CTD  stations 
occupied  in  the  northern  and  central  parts  of  the  basin.  As  it  can  be 
seen  from  Fig.  1,  most  of  the  CTD  soundings  (216)  were  shallower 
than  280  m,  with  only  a  three  casts  taken  at  deeper  locations.  The 
total  number  of  TS  observations  used  in  the  assimilation  is  9650.  With 
the  total  number  of  the  observed  velocities  13,856  the  dimension  of 
the  observation  space  was  I<  =  23,506. 

2.3.  Assimilation  parameters 

In  the  course  of  the  experiments  described  in  the  next  section,  we 
tried  to  keep  the  parameters  of  the  tested  4dVar  and  a4dVar  systems 
as  close  as  possible  to  each  other.  However,  due  to  the  different  for¬ 
mulations  (observation  space  for  4dVar  and  state  space  for  a4dVar), 
certain  discrepancies  remained  in  the  shape  of  the  background  error 
covariance  B.  In  both  algorithms  B  is  represented  by  the  product  VCV 
where  V  is  the  diagonal  matrix  of  the  background  error  rms  variances 
and  C  is  the  respective  correlation  matrix. 

In  the  4dVar  algorithm,  C  is  represented  implicitly  by  the  kernel 
of  the  heat  transfer  equation 

C~exp^a2A^,  (8) 

where  a  is  the  decorrelation  length  scale  and  A  is  the  discretized  2d 
Laplacian  operator.  Numerically,  the  action  of  C  on  a  state  vector  is 
computed  by  integrating  the  heat  transfer  equation  (e.g.,  Weaver  and 
Courtier,  2001 ).  In  the  vertical,  the  decorrelation  scale  was  set  to  zero. 
The  correlation  matrix  (8)  is  rank-deficient,  so  the  4dVar  solution  is 
obtained  in  the  range  of  B. 
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Fig.  3.  The  background  error  correlation  functions. 


The  a4dVar  algorithm  is  formulated  in  the  state  space  and  requires 
a  definition  of  C-1  which  was  explicitly  specified  as  the  inverse  of  the 
second-order  approximation  of  the  exponent  in  (8) 


C 


-l 

a 


(9) 


where  I  is  the  identity  operator  in  state  space,  and  b  =  8/Tta  to  pre¬ 
serve  the  value  of  the  integral  decorrelation  scale  specified  in  4dVar 
(e.g.,  Yaremchuk  and  Smith,  2011). 

Although  the  respective  correlation  functions  are  somewhat  dif¬ 
ferent  in  shape  (Fig.  3),  we  assume  that  this  difference  has  minor  ef¬ 
fects  on  the  overall  results  of  the  assimilation  compared  to  the  effect 
of  non-linearity. 

The  rest  of  the  assimilation  parameters  were  identical  for  both  the 
4dVar  and  a4dVar  assimilation  systems.  The  value  of  a  was  chosen  to 
be  9  km,  consistent  with  typical  estimates  of  the  baroclinic  deforma¬ 
tion  radius  in  the  region  (e.g.,  Cushman-Roisin  and  Korotenko,  2007). 
The  background  error  rms  variances  (diagonal  elements  of  V)  were 
assumed  to  be  proportional  to  the  rms  variability  of  the  respective 
NCOM  fields  from  the  first  guess  solution  with  the  typical  errors  of 
1  °C,  0.1  psu,  and  10  cm/s  near  the  surface.  Observation  errors  were 
assumed  to  be  spatially  uncorrelated  with  the  rms  variances  (diago¬ 
nal  elements  of  Rl/2)  dependent  only  on  depth  in  the  manner  shown 
by  the  solid  profiles  in  Fig.  2).  Actual  vertical  distributions  of  the  ob¬ 
servational  rms  error  variances  were  specified  by  multiplying  these 
curves  by  0.33  for  temperature  and  salinity,  and  by  0.5  for  velocity. 

The  a4dVar  EOF  analysis  was  performed  with  respect  to  the  diag¬ 
onal  metric  specified  by  the  inverse  background  error  variances  V  2. 

The  stopping  criteria  for  the  iterative  processes  were  selected  as 
follows:  for  the  4dVar  system  the  solution  of  the  system  for  the  rep¬ 
resenter  coefficients  was  terminated  after  nt=  7  iterations,  when  the 
accuracy  of  the  conjugate  gradient  (CG)  solver  was,  as  a  rule,  below 
10~3.  Experiments  with  the  larger  number  of  CG  iterations  (inner 
loops)  have  shown  only  minor  effects  on  the  final  optimal  solution, 
whereas  reduction  of  nt  resulted  in  amplification  of  the  effects  related 
to  the  instabilities  of  the  tangent  linear  model  and/or  its  adjoint,  caus¬ 
ing  a  sharp  increase  of  the  cost  function  after  3-4  outer  loops.  With 
the  value  of  nt= 7,  8-10  outer  loops  were  executed  before  the  values 
of  J  started  to  increase. 

For  the  a4dVar  system,  the  minimization  was  terminated  when 
the  total  CPU  time  reached  the  value  used  by  the  respective  4dVar 
experiment.  The  number  of  ensemble  members  was  kept  constant  at 
ms  =  9  through  all  the  experiments. 

An  important  technical  issue  in  a4dVar  was  the  choice  of  the  per¬ 
turbation  magnitude  e.  Ideally,  this  value  should  be  as  small  as  pos¬ 
sible  to  keep  the  perturbed  system  linear.  In  practice,  values  of  e  sig¬ 
nificantly  below  10-2  were  ineffective  due  to  the  loss  of  accuracy  in 
calculating  the  perturbations  of  Y,  ,  especially  in  the  course  of  14  days 
of  model  integration.  For  that  reason,  on  each  iteration,  the  value  of 
e  was  selected  in  a  way  that  only  one  field  among  all  the  ms  pertur¬ 
bations  pj  could  reach  its  critical  magnitude  at  one  point  of  the  do¬ 
main.  The  respective  critical  values  of  the  temperature,  salinity,  and 


velocity  perturbations  were  set  to  4  °C,  2  psu,  and  0.5  m/s,  respec¬ 
tively.  This  strategy  allowed  fast  convergence  while  avoiding  devel¬ 
opment  of  instabilities  within  the  perturbed  model  runs. 

Preliminary  experiments  were  also  performed  to  tune  the  number 
of  search  subspaces  k  to  be  kept  in  memory  along  with  the  vectors  Y, 
needed  for  H-orthogonalization.  The  large  number  of  elements  in  Y 
made  the  orthogonalization  process  rather  time-consuming  for  k  > 
10.  Besides,  H-orthogonality  was  quickly  lost  with  iterations  due  to 
strong  non-linearity  of  the  model  and  the  above  mentioned  inaccu¬ 
racy  in  estimating  the  perturbations  of  Y  due  to  the  finite  value  of  e. 
After  some  experimentation,  it  was  found  that  k  =  3  with  the  above 
mentioned  strategy  of  selecting  e  provided  the  fastest  convergence 
for  the  a4dVar  algorithm. 

The  EOF  reduction  of  the  search  subspace  (step  2  in  the  lay¬ 
out  of  Section  2.1.2)  may  seem  to  be  redundant  for  a  linear  system, 
but  appears  to  be  important  in  the  considered  a4dVar  application: 
First,  performing  the  search  along  the  few  principal  modes  extracted 
from  the  time  sequence  {rn}  tends  to  keep  the  minimization  process 
within  the  most  persistent  (geostrophically  and  hydrostatically  bal¬ 
anced)  manifold,  thus  avoiding  searches  over  initial  states  that  tend 
to  generate  excessive  ageostrophic  (i.e.,smaller-scale)  motions.  Sec¬ 
ond,  rescaling  the  EOF  metric  V  proved  to  be  useful  in  restarting  the 
H-orthogonalization:  rescaling  was  done  every  time  when  the  rela¬ 
tive  reduction  y  =  SJ/J  of  the  cost  function  at  the  start  of  the  new 
orthogonalization  cycle  was  ten  times  smaller  than  the  mean  value 
of  y  on  the  previous  cycle.  In  the  event  y  <  0.1  y,  salinity  entries  of 
V  were  inflated  by  the  factor  of  5  and  then  restored  to  their  original 
values  on  the  next  occurrence  of  the  event. 

3.  Results 

In  the  reported  experiments  we  varied  the  length  of  the  assimila¬ 
tion  window  from  short  (4  days,  N  =  9)  to  moderate  (8  days,  N  =  17) 
and  long  (14  days,  N= 29)  duration.  Performance  of  the  assimilation 
algorithms  was  evaluated  in  three  categories:  the  forecast  skill  at  the 
end  of  the  assimilation  window  (for  N  =  917),  the  rate  of  conver¬ 
gence,  and  by  qualitative  inspection  of  the  optimal  model  trajectories. 

3.1.  Convergence  rates  and  computational  expense 

To  assess  the  rates  of  convergence,  one  has  to  have  an  ability  to 
compare  the  reduction  of  the  cost  function  with  iterations,  which  is 
not  straightforward  for  two  reasons. 

First,  in  the  4dVar  algorithm  considered  here,  the  regularization 
term  of  the  cost  function  can  be  evaluated  only  within  the  range  of 
the  correlation  matrix  defined  by  (9).  To  avoid  the  burden  of  restrict¬ 
ing  CQ  to  the  range  of  C,  we  compared  only  the  observational  parts  of 
the  4dVar  and  a4dVar  cost  functions  (second  term  in  Eq.  (3)). 

Second,  the  number  of  iterations  required  for  convergence  cannot 
be  considered  as  an  objective  criterion  because  4dVar  and  a4dVar  it¬ 
erations  are  different  in  nature.  Due  to  the  non-linearity  of  the  prob¬ 
lem,  an  iteration  (either  4dVar  or  a4dVar)  performs  minimization  in 
the  vicinity  of  the  current  (suboptimal)  state,  but  4dVar  does  that  in 
the  range  of  B,  whereas  a4dVar  minimizes  in  the  subspace  of  a  much 
smaller  dimension  spanned  by  pm.  For  that  reason,  iterations  require 
quite  different  computational  resources  and  should  be  compared  in 
terms  of  CPU  time. 

Fig.  4  shows  such  a  comparison  by  rescaling  the  horizontal  axis 
with  the  total  CPU  time  za  required  by  one  a4dVar  iteration.  The 
value  of  ra  was  11  times  larger  than  the  CPU  time  rm  of  a  direct 
NCOM  model  run  for  a  given  experiment,  i.e.  ta  ~  11  rm.  The  major 
contribution  to  ra  is  given  by  the  ensemble  run  (9rm,  p.3  in  the  lay¬ 
out  of  Section  2.1.2),  while  the  master  NCOM  run  (p.l )  and  operations 
listed  in  pp.2  and  4  require  rm  and  0.8rm,  respectively.  Overall,  con¬ 
vergence  was  achieved  at  an  expense  of  60-70  iterations  (650-800 
NCOM  runs). 
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Fig.  4.  Relative  reduction  i]  of  the  cost  function  with  iterations  (marked  by  circles)  for 
different  assimilation  periods.  The  horizontal  axis  is  scaled  by  the  CPU  time  required 
for  the  a4dVar  iteration.  Squares  label  the  4dVar  outer  loops. 


As  may  be  seen  in  Fig.  4,  a  single  4dVar  iteration  was  approxi¬ 
mately  equivalent  to  6-7  a4dVar  iterations,  or  70-80  direct  model 
runs.  This  computational  expense  arises  because  sequential  execu¬ 
tion  of  the  adjoint  and  tangent  linear  codes  (inner  loops  of  the  CG 
solver)  required  around  llrm,  whereas  one  4dVar  outer  loop  in¬ 
cluded  seven  inner  loops  to  solve  the  system  of  linear  equations  for 
the  representer  coefficients.  In  a  series  of  experiments,  it  was  found 
that  executing  seven  inner  loops  provided  a  103-fold  reduction  of  the 
system’s  residual  and  was  optimal  with  regard  to  the  total  CPU  time 
of  the  4dVar  algorithm. 

Fig.  4  shows  that,  in  general,  the  tested  a4dVar  method  is  com¬ 
putationally  comparable  to  the  observation  space  4dVar.  Although 
the  total  CPU  time  required  for  reduction  of  J  by  the  factor  of  0.4 


(attained  after  the  first  outer  loop  of  the  4dVar)  appears  to  be  similar 
for  the  4dVar  and  a4dVar  methods,  the  a4dVar  minimization  notice¬ 
ably  slows  down  at  subsequent  iterations,  especially  for  longer  assim¬ 
ilation  windows  (r=8,14  days).  This  could  be  partly  explained  by  the 
fact  that,  with  longer  windows,  operators  M"  tend  to  depart  farther 
away  from  the  identity  and  it  becomes  increasingly  more  difficult  for 
the  a4dVar  algorithm  to  find  the  minimum  without  the  additional 
information  on  the  structure  of  H  provided  by  the  adjoint  code  in 
4dVar. 

Fig.  5  demonstrates  the  time  evolution  of  the  quantities 

/,"  =  ([(Hnx,n  -  d^)T(Hnx2  -  dq")/nq]1/2J  (10) 

characterizing  the  daily  averaged  ()  model-data  misfits  of  the  vari¬ 
ous  state  vector  components  before  (black  lines)  and  after  (gray  lines) 
optimization  with  a  14-day  assimilation  window  (i.e.,  using  all  the 
available  data).  The  subscript  q  takes  the  values  of  the  labels  in  the 
mid-bottom  parts  of  Fig.  5  which  indicate  the  observed  variables 
(temperature,  salinity  and  velocity  vector)  for  which  the  statistics  fq 
were  computed,  whereas  nq  stands  for  the  total  number  of  respective 
observations  taken  at  a  given  day. 

Comparison  of  the  model-data  differences  for  the  background 
forecast  (thick  black  lines  in  Fig.  5)  with  those  computed  for  persis¬ 
tence  (thin  lines)  shows  their  approximate  similarity,  especially  for 
fT  and  fs.  This  similarity  can  be  partly  explained  by  small  biases  in 
the  temperature  and  salinity  fields  of  the  background  solution  and 
large  discrepancies  in  representation  of  the  mesoscale  structures  by 
the  background  solution.  As  a  consequence,  persistence  assumption 
appears  to  be  much  less  valid  for  the  velocity  field  (Fig.  5,  lower  panel) 
which  provides  the  major  contribution  to  the  combined  behavior  of 
fq  shown  in  the  upper  panel. 

The  upper  panel  in  Fig.  5  also  shows  a  remarkable  similarity  in  the 
time  evolution  of  the  combined  model-data  misfit  for  the  4dVar-  and 
a4dVar-optimized  NCOM  states.  The  a4dVar  algorithm  has,  however, 
a  noticeable  tendency  to  provide  a  better  fit  at  the  beginning  of  the 
assimilation  window,  clearly  visible  in  the  lower  panels  for  fs  and  fv. 
This  can  be  explained  by  the  above  mentioned  property  of  a4dVar  to 
better  retrieve  optimal  states  at  shorter  integration  times. 

When  separated  into  different  components,  behavior  of  ff,  /", 
and  fj]  reveals  more  differences.  In  particular,  the  4dVar  method  pro¬ 
vides  a  much  better  fit  to  the  temperature  data  after  August  20  (in 
the  second  half  of  the  assimilation  window),  but  appears  to  be  10- 
13%  worse  than  a4dVar  with  respect  to  salinity  and  velocity  data. 

A  large  contribution  to  a  better  salinity  fit  is  given  by  the  first  two 
days  of  the  a4dVar  model  trajectory  (third  panel  in  Fig.  5).  However, 
certain  gains  relative  to  4dVar  are  also  observed  at  the  end  of  the 
assimilation,  which  is  quite  opposite  to  the  difference  in  the  values  of 
fT.  We  attribute  the  better  salinity  fit  to  the  variable  EOF  metric  used 
during  generation  of  the  a4dVar  search  directions  (Section  2.3). 

Compared  to  fT  and  fs,  the  overall  improvement  of  the  model-data 
misfit  is  the  smallest  for  velocity  (bottom  panel  in  Fig.  5),  which  was 
characterized  by  the  observation  errors  of  R1/2  ~  7-10  cm/s  in  the 
cost  function.  Several  assimilation  runs  were  made  with  significantly 
smaller  (3-5  cm/s)  errors,  but  they  were  found  to  be  inconsistent 
with  a  posteriori  statistics  of  the  model-data  misfits  as  the  optimal 
cost  function  values  in  these  cases  were  much  larger  than  those  ob¬ 
tained  in  the  reported  experiments.  The  a4dVar-optimized  value  of  fv 
is  persistently  smaller  during  the  entire  assimilation  period  providing 
the  13%  better  value  (as  compared  to  4dVar)  in  the  14-day  average. 
This  advantage  could  be  partly  attributed  to  the  fact  that  the  a4dVar 
search  directions  are  derived  from  the  most  persistent  patterns  of 
the  model-data  misfits  and  therefore  tend  to  be  closer  to  the  slowly 
evolving  (geostrophically  and  hydrostatically  balanced)  modes  of  the 
flow.  This  property  was  also  reflected  in  the  better  velocity  forecast 
skill  of  the  a4dVar  solutions. 
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Fig.  5.  Evolution  of  the  root-mean-square  model-data  misfits  /,  characterizing  the 
background  (BG,  thick  black  lines),  4dVar-optimized  (4d,  thin  dashed  gray  lines)  and 
a4dVar-optimized  (a4d)  solutions.  Thick  dash-dotted  line  shows  the  misfit  with  the 
background  fields  at  t  =  0  (persistence).  The  values  of/,  are  shown  on  the  right  axis  of 
each  panel.  The  left  axis  quantifies  the  number  of  the  data  points  for  each  day  in  thou¬ 
sands  (shown  by  gray  shaded  rectangles).  The  ratio  of  the  mean  values  of  /,  averaged 
over  the  assimilation  window  for  the  4dVar  and  a4dVar  methods  is  given. 

3.2.  Forecast  skill 

The  quality  of  the  assimilated  solutions  was  assessed  for  4-  and 
8-day  experiments  using  comparison  with  observations  outside  the 
respective  assimilation  windows.  Evolution  of  the  quantities/1  for  the 
background  and  optimized  solutions  is  shown  in  Fig.  6  for  the  4-day 
assimilation  experiment. 

The  general  behavior  of /is  consistent  with  the  one  obtained  in  the 
14-day  experiment,  showing  persistently  better  4dVar  forecasts  in 
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Fig.  6.  Forecast  skills /s,  /„  and/t  of  the  4dVar  and  a4dVar-optimized  solutions  for  the 
4-day  assimilation  window.  The  relative  number  of  the  respective  data  points  for  each 
day  is  shown  by  gray  shaded  rectangles.  Vertical  dashed  line  show  the  time  interval  of 
data  assimilation.  Ratios  of  the  mean/values  averaged  over  the  3-  and  9-days  intervals 
are  shown. 


temperature  and  the  advantage  of  a4dVar  in  the  salinity  and  velocity 
forecasts.  The  upper  panel  in  Fig.  6  summarizes  the  forecast  skill  and 
indicates  that  4dVar  slightly  (4-5%)  outperforms  a4dVar,  mostly  be¬ 
cause  of  the  better  temperature  forecasts  (second  panel  from  above). 
On  the  other  hand,  the  4dVar-optimized  salinity  is  characterized  by 
very  low  forecast  skill,  especially  during  August  21-25,  when  it  was 
even  farther  away  from  the  observations  than  the  background  fore¬ 
cast. 

The  4dVar-optimized  velocities  show  only  small  improvements 
compared  to  the  background  solution  (bottom  panel  in  Fig.  6).  In  con¬ 
trast,  the  a4dVar-optimized  velocities  demonstrate  10-30%  reduction 
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of  the  model-data  misfit  within  the  assimilation  window,  which  per¬ 
sists  for  up  to  three  days  (August  18-21 )  of  the  free  model  run.  Af¬ 
ter  August  21,  the  velocity  mismatch  of  the  background,  a4dVar  and 
4dVar-optimized  solutions  are  nearly  identical.  Qualitatively  similar 
behavior  of  the  forecast  skill  and  its  distribution  among  the  state  vec¬ 
tor  components  was  observed  in  the  results  of  the  8-day  assimilation 
experiment. 

In  interpreting  the  forecast  skill  assessment,  it  is  necessary  to  take 
into  account  that  the  temperature  and  salinity  observations  after  Au¬ 
gust  21  are  much  less  numerous  than  those  taken  between  the  18 
and  21  of  August  (cf.  gray  rectangles  of  the  second  and  third  panels  in 
Fig.  5),  whereas  the  velocity  data  cover  a  rather  limited  area  shown 
by  the  triangles  in  Fig.  1. 

In  general,  the  overall  forecast  skill  provided  by  the  a4dVar 
method  appears  to  be  comparable  with  that  of  the  4dVar  (upper 
panel  in  Fig.  6),  and  in  some  aspects  (such  as  short-term  velocity 
forecast),  the  a4dVar  technique  provides  noticeably  better  results. 
It  should  be  noted  that  available  observations  could  effectively  con¬ 
strain  only  a  small  part  K/M  =23,506/1,493,570  ~  1.5%  of  the  model’s 
degrees  of  freedom,  so  one  should  expect  substantial  differences  in 
the  small-scale  structure  of  the  optimal  solutions  obtained  by  two 
different  methods. 

3.3.  Comparison  of  the  optimal  solutions 

Temperature  and  velocity  increments  for  the  optimal  states  of  the 
14-day  assimilation  experiment  are  shown  in  Fig.  7.  A  certain  coher¬ 
ence  between  the  larger  scale  corrections  to  the  background  temper¬ 
ature  field  are  clearly  seen  in  the  northern  part  of  the  model  domain 
that  is  well  covered  by  observations  (cf.  Fig.  1 ).  The  time-mean  corre¬ 
lation  coefficients  p  between  the  low-pass  filtered  temperature  and 
salinity  increments  of  the  4dVar  and  a4dVar  solutions  are  0.61  and 
0.45,  respectively  if  averaging  is  performed  in  the  upper  200  m  over 
the  northern  part  of  the  domain.  In  the  data-free  region  south  of  the 
340  km  mark,  the  correlations  are  substantially  lower  (respectively, 
0.26  and  0.32)  and  lie  below  the  95%  confidence  level  of  nonzero  cor¬ 
relation  (0.36).  Similar  values  of  p  (0.59  and  0.32  in  the  northern  and 
southern  subregions,  respectively)  were  obtained  for  the  sea  surface 
height  field. 

Velocity  increments  appear  to  have  the  lowest  correlations  among 
the  model  fields  with  time-averaged  values  of  pv  =0.36,  0.27  for  the 
northern  and  southern  subregions,  respectively.  The  lowest  corre¬ 
lations  [pv  =  0.09,  pT  =  0.21,  and  ps  =  0.12)  were  observed  in  the 
data-free  southern  subregion  during  the  first  4  days  (8/14-8/18  )  of 
the  assimilation.  Such  incoherence  between  the  increments  is  caused 
by  excessive  ageostrophic  activity  (left  panel  in  Fig.  7)  of  the  4dVar 
solution  at  the  beginning  of  the  assimilation  window.  Ageostrophic 
nature  of  this  disturbances  is  quantified  by  4  times  larger  magni¬ 
tude  of  the  divergence  field  as  compared  to  the  rest  of  the  domain. 
The  ageostrophic  mode  disappears  at  the  later  times  and  does  not 
affect  the  cost  function  because  the  southern  subregion  is  virtually 
data-free,  whereas  smoothness  constraints  are  imposed  on  the  model 
fields  only  at  the  initial  time. 

The  problem  could  be  solved,  apparently,  by  introducing  balance 
constraints  (e.g.,  Weaver  et  al.,  2005)  into  the  definition  of  the  back¬ 
ground  error  covariance  at  n  =  0,  which  may  not  be  necessary  if  the 
NCOM  4dVar  were  run  in  the  weakly  constrained  mode,  i.e.,  if  back¬ 
ground  errors  were  prescribed  throughout  the  entire  assimilation 
window.  For  comparison  purposes  we  ran  the  4dVar  system  in  the 
strongly  constrained  mode  and  the  effect  became  visible  after  several 
outer  loops.  It  is  remarkable  that  the  a4dVar  algorithm  appears  to  be 
much  less  susceptible  to  excitation  of  the  ageostrophic  modes  (right 
panel  in  Fig.  7),  possibly  because  the  EOF-derived  descent  directions 
span  subspaces  characterized  by  slower  time  variation  of  the  model 
trajectory  and,  therefore,  tend  to  be  closer  to  geostrophic  and  hydro¬ 
static  balance.  The  null-space  nature  of  the  ageostrophic  features  at 


Fig.  7.  Temperature  and  velocity  differences  between  the  background  and  optimized 
NCOM  states  at  20  m  on  August  15  03  UTC.  Results  of  4dVar  and  a4dVar  optimiza¬ 
tions  are  shown  in  the  left  and  right  panels,  respectively.  Velocity  scale  is  shown  at  the 
bottom  of  the  left  panel. 


Fig.  8.  Distances  between  the  background  (BG)-a4dVar  (thick  black  lines),  BG-4dVar 
(gray  lines),  and  a4dVar-4dVar  (thin  black  lines)  solutions  after  the  first  4dVar  outer 
loop  (1),  third  loop  (3)  and  upon  convergence  (oo).  Respective  4dVar  outer  loops  are 
labeled  by  squares  in  the  middle  and  lower  panels  of  Fig.  4.  The  a4dVar  distances  are 
given  on  the  iterations  (circles  in  Fig.  4)  corresponding  to  the  equivalent  CPU  times  of 
the  respective  4dVar  outer  loops  (squares  in  Fig.  4).  All  the  distances  are  normalized 
by  the  magnitude  of  the  background  trajectory.  Angles  are  assessed  using  the  scalar 
product,  associated  with  Eq.  (12). 


the  beginning  of  the  4dVar-optimized  model  trajectory  can  be  also 
traced  by  looking  at  the  evolution  with  iterations  of  the  distances  Sx 
between  the  background  and  (sub)optimal  model  trajectories  (Fig.  8): 
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As  can  be  seen,  the  largest  (60-70%)  a4dVar  correction  of  the  back¬ 
ground  trajectory  occurs  at  the  first  few  iterations,  whereas  4dVar  at 
the  first  iteration  makes  a  significantly  smaller  correction  and  subse¬ 
quently  produces  larger  updates  that  are  characterized  by  a  relatively 
small  reduction  of  the  cost  function  (cf.  Fig.  4).  Such  a  behavior  is 
typical  for  a  null  space  search  which  appears  to  be  inhibited  in  the 
a4dVar  case.  Introduction  of  the  balance  constraints  into  B  will  cer¬ 
tainly  improve  the  performance  of  both  algorithms  with  a  potentially 
larger  benefit  for  the  4dVar  case. 

4.  Discussion  and  conclusions 


The  major  goal  of  the  present  study  was  to  show  the  feasibility 
of  the  a4dVar  technique  (Yaremchuk  et  al.,  2009)  in  a  realistic  ap¬ 
plication  and  compare  its  performance  with  the  4dVar  method.  We 
have  found  that  the  adjoint-free  approach  is  capable  of  producing  op¬ 
timized  solutions  of  similar  quality  to  4dVar  with  comparable  com¬ 
putational  expense.  It  was  also  found  that  the  a4dVar  technique  is 
less  susceptible  to  excitation  of  ageostrophic  modes  in  the  data-free 
regions  if  balance  constraints  are  not  imposed  on  the  background  er¬ 
ror  covariances. 

A  distinctive  feature  of  the  a4dVar  technique  presented  here  is  the 
iterative  approach  to  minimization  of  the  cost  function.  The  adjoint- 
free  methods  of  Qui  et  al.  (2007)  and  Liu  et  al.  (2008)  employ  mini¬ 
mization  in  a  predetermined  subspace  derived  from  the  background 
error  statistics  (with  or  without  localization  of  the  background  er¬ 
ror  covariance).  The  tested  a4dVar  method  follows  the  same  princi¬ 
ple  as  the  model-data  misfits  are  projected  on  the  range  of  B  (Eq.  (7)), 
nonetheless  convergence  to  xopt  cannot  generally  be  guaranteed  since 
the  process  may  be  subject  to  breakdown  (situations  when  new 
search  directions  are  linearly  dependent  on  previous  ones).  In  future 
applications,  the  issue  could  possibly  be  resolved  using  the  meth¬ 
ods  already  developed  for  the  GMRES-type  algorithms  (e.g.,  Reichel 
and  Ye,  2005),  which  are  not  immune  to  breakdowns,  similar  to  the 
a4dVar  minimization  technique.  Experiments  reported  here  do  show, 
however,  that  the  observed  rate  of  convergence  of  the  a4dVar  mini¬ 
mization  process  may  slow  down  relative  to  what  is  seen  with  4dVar 
and  this  effect  may  be  exacerbated  with  increasing  size  of  assimi¬ 
lation  window  (Fig.  4).Nonetheless,  a4dVar  tends  to  produce  better 
results  at  the  earlier  stages  of  assimilation  than  4dVar  and  in  general 
its  performance  could  still  be  viewed  as  satisfactory. 

Taking  advantage  of  the  trend  toward  massive  parallelization  in 
computer  technologies,  the  adjoint-free  variational  methods  esti¬ 
mate  the  cost  function  gradient  with  a  “brute  force"  approach  that 
employs  finite  difference  approximations  along  predetermined  direc¬ 
tions  in  state  space.  Selection  of  the  most  effective  directions  (search 
subspaces  Sf)  based  on  this  sampling  becomes  an  issue  of  primary 
importance.  Experience  shows  that  there  exists  a  considerable  free¬ 
dom  in  generating  search  directions  (Yaremchuk  et  al.,  2009;  Pan¬ 
teleev  et  al„  2015)  as  long  as  they  are  kept  being  spatially  smooth  and 
H-orthogonal.  In  particular,  building  §,  on  the  eigenvectors  of  B  in 
the  descending  order  of  their  eigenvalues  may  also  work  reasonably 
well,  as  shown  in  Appendix  A.  In  the  reported  experiments,  we  inves¬ 
tigated  a  number  of  methodologies  in  building  §,  and  found  Eq.  (7)  to 
be  the  most  effective  computationally.  This  methodology  can  be  de¬ 
veloped  further  by  introducing  balance  constraints  into  B_1.  In  this 
case  the  inverse  background  error  covariance  should  be  replaced  by 
the  composite  matrix 
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where  B^1  and  B71  are  the  inverse  covariances  of  the  unbalanced 
components  of  the  state  vector  (e.g.,  defined  by  (9))  and  L  is  the  bal¬ 
ance  operator.  Further  improvements  can  be  made  by  replacing  the 
diffusion  operator  in  Eq.  (9)  with  a  more  general  expression  (e.g., 
Weaver  and  Mirouze,  2012;  Yaremchuk  and  Nechaev,  2013). 


An  important  issue  with  the  a4dVar  technique  is  its  extension  to 
optimization  of  other  sets  of  variables  that  may  control  the  model 
trajectory,  such  as  surface  forcing  fields.  One  of  the  possible  solu¬ 
tions  in  this  case  augments  the  search  subspaces  (ocean  model  states) 
by  the  leading  EOFs  of  the  surface  forcing  error  fields.  This  will  re¬ 
quire  a  better  knowledge  of  error  statistics  of  the  atmospheric  model 
used  to  force  the  ocean  in  a  particular  application.  In  view  of  recent 
rapid  development  of  the  observational  systems  and  data  acquizition 
techniques  in  the  atmosphere,  the  issue  of  accessibility  to  the  above 
mentioned  statistics  seems  likely  to  be  resolvable  in  the  near  term. 
Moreover,  the  a4dVar  technique  appears  to  be  even  more  suitable 
for  coupled  ocean-atmosphere  systems,  where  external  forcing  er¬ 
rors  tend  to  play  a  lesser  role  at  the  time  scale  of  a  typical  assimilation 
window. 

In  terms  of  the  computational  expense,  the  a4dVar  technique  ap¬ 
pears  roughly  comparable  to  4dVar,  mostly  because  of  the  excessive 
computational  cost  of  tangent  linear  and  adjoint  codes  that  were,  on 
average,  several  times  more  expensive  than  a  direct  run  of  the  non¬ 
linear  NCOM  model  (a  typical  situation  with  state-of-the-art  OGCMs, 
e.g.,  Oldenborgh  et  al.,  1999  and  Heimbach  et  al.,  2005).  On  massively 
parallel  machines,  the  advantage  of  a4dVar  will  be  more  noticeable 
due  to  the  limited  parallel  scalability  of  an  OGCM  code,  be  it  original, 
adjoint,  or  tangent  linear. 

A  much  larger  computational  advantage  is  evident  when  consider¬ 
ing  the  wall  time  in  a  massively  parallel  environment,  which  formally 
allows  a4dVar  to  search  over  multiple  directions  at  a  fraction  of  the 
wall  time  used  by  4dVar  to  generate  a  steepest  descent  direction.  In 
fact,  in  the  reported  experiments  with  NCOM  model,  one  a4dVar  run 
was  executed  almost  five  times  faster  if  all  the  ensemble  members 
were  run  on  separate  nodes.  This  property  of  the  a4dVar  approach 
gives  good  prospects  for  its  further  development  in  sync  with  and 
other  types  of  ensemble  data  assimilation  techniques  that  are  based 
on  relaxed  communication  requirements  between  processors. 
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Appendix  A.  4dVar/a4dVar  comparison  in  the  linear  case 

To  compare  the  performance  of  4dVar  and  a4dVar  techniques  in 
a  simple  environment,  consider  the  problem  of  retrieving  the  initial 
field  of  tracer  concentration  rj(x,  0)  from  observations  at  some  dis¬ 
tant  time  T.  The  tracer  evolution  is  governed  by 

dtil  +  uV?7  -  fiAr]  =  f(x,  t)  (A.l) 

in  a  closed  rectangular  49  x  91  domain  £2  with  the  boundary  con¬ 
dition  ii(dS2,  t)  =  0.  Eq.  (A.l)  is  discretized  on  a  regular  grid  using 
first-order  RK  time-stepping,  upwind  advection,  and  a  standard  5- 
point  stencil  for  the  Laplacian  with  unit  steps  in  temporal  and  spatial 
directions.  Velocity  u  =  (u,  v)  at  any  space-time  location  was  defined 
by  u  =  -0.2  +  0.01  u;  v  =  -0.1  +  0.01  v,  where  v  is  the  white  noise  on 
unit  interval.  The  forcing/was  generated  by  setting  / (x,  t)  =  .001  v  in 
every  point  of  the  space-time  grid.  The  coefficient  fi  was  set  to  10~5, 
so  that  diffusion  was  largely  determined  by  the  numerics. 

The  simulated  data  comparison  experiment  was  set  as  follows. 
Given  the  initial  tracer  distribution  fj  =  r){x,  0)  =  exp[-(x  -  x0)2/9] 
with  x0  =  (70,35)  (bell-shaped  disturbance  in  Fig.  Ala),  the  model 
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Fig.  Al.  Reconstruction  of  the  initial  condition  of  the  tracer  field  by  4dVar  (b)  and 
a4dVar  (c,d)  techniques.  Composite  map  of  the  reconstructed  tracer  field  evolution  is 
shown  in  the  upper  panel,  (a):  Initial  position  of  the  reconstructed  feature  (Gaussian 
eddy  at  x  =  70,  y  =  35  km)  is  superimposed  on  the  tracer  field  (contours)  at  the  ob¬ 
servation  time  t  =  200  when  the  eddy  diplaced  to  x  =  25,  y  =  15  km.  Circles  denote 
observation  points.  The  errors  in  approximation  of  the  true  field  at  t  =  0  are  shown  in 
the  left  corners  of  the  panels  b-d,  showing  the  initial  field,  reconstructed  by  various 
methods. 


was  integrated  for  T  =  200  time  steps  to  obtain  the  final  distribu¬ 
tion  pix.T)  shown  by  contours  in  the  same  panel.  The  initial  dis¬ 
turbance  almost  completely  dispersed  and  migrated  to  xT  ~  (25, 15) 
(see  contours  in  the  same  panel).  After  that,  rj(x,  T)  was  sampled  at 
200  points  shown  in  Fig.  Ala,  and  obtained  numbers  were  used  to  re¬ 
construct  fj  by  minimizing  the  cost  function  (3)  under  the  dynamical 
constraint  (A.l )  with  the  background  error  covariance  defined  by  (9). 

Optimal  approximations  fj-fj  obtained  by  4dVar  and  a4dVar  tech¬ 
niques  are  shown  in  Fig.  Alb  and  Fig.  Ale,  respectively).  To  specify 
search  directions  in  the  a4dVar  method,  200  observations  were  split 
into  ms  =  10  equal  groups  so  an  observation  operator  for  each  search 
direction  in  (7)  had  20  observation  locations. 

The  quality  of  reconstruction  was  assessed  by  the  parameter 

e  =  V((fj-ri)2)/(f}2)  (A.2) 

where  angular  brackets  denotes  averaging  over  small  rectangles  in 
Fig.  Al.  Comparison  of  Fig.  Alb  and  c  suggests  that  the  a4dVar 
method  is  capable  of  providing  a  solution  of  the  same  quality  with 
4dVar. 

Fig.  Aid  illustrates  another  a4dVar  solution,  using  search  sub¬ 
spaces  specified  by  the  eigenvectors  of  B  in  the  decreasing  order  of 
their  eigenvalues  (in  this  case,  sines  with  decreasing  wavelengths  in 
both  directions).  The  result  obtained  is  of  similar  quality,  suggesting 
that  the  general  a4dVar  strategy  of  minimizing  J  over  a  sequence  of 


Fig.  A2.  Reduction  of  the  cost  function  against  CPU  time  for  4dVar  and  a4dVar  tech¬ 
niques.  The  4dVar  CPU  time  is  multiplied  by  five  to  mimic  larger  CPU  requirements 
of  the  state-of-the-art  adjoint  models.  Inset:  Convergence  of  the  a4dVar-B  solution 
(Fig.  Aid)  to  the  exact  solution.  Dashed  line  shows  the  convergence  rate  given  by  (A.5). 


smooth  H-orthogonal  SDs  may  work  well  with  various  methods  of 
generating  search  directions. 

In  terms  of  computational  expense,  the  4dVar  method  pro¬ 
vided  approximately  five  times  faster  reduction  of  the  cost  function 
(Fig.  A2)  due  to  high  efficiency  of  the  adjoint  model.  In  this  simple 
case,  an  adjoint  model  run  required  the  same  amount  of  time  as  the 
direct  model  run.  In  real  applications,  the  tangent  linear  and  adjoint 
codes  are  several  times  more  expensive  to  run  and  the  a4dVar  tech¬ 
niques  may  prove  to  be  more  competitive,  as  shown  in  Section  3  of 
the  present  study. 

Finally,  the  known  spectrum  of  B  provides  an  opportunity  to  as¬ 
sess  the  convergence  rate  of  the  ad4Var  solution  exposed  in  Fig.  Aid. 
Assume  that  after  k  a4dVar  iterations  m  =  kms  H-orthogonal  direc¬ 
tions  have  been  already  searched  and  the  kth  approximation  fjk  to 
the  optimal  solution  fj  =  H  1  b  have  been  found.  Without  loss  of 
generality,  the  eigenvectors  0,  of  B  could  be  normalized  to  satisfy 
0jB~'0j  =  1,  so  that  their  (Euclidean)  norm  is  equal  to  the  associ¬ 
ated  eigenvalue  cr,.  The  magnitude  em  of  the  approximation  error 
em  =  fj  -  fjm  with  respect  to  the  norm  induced  by  the  inverse  covari¬ 
ance  can  be  assessed  by  projecting  fj  on  the  unexplored  directions: 

OO 

em=e^B-,em<  |t7TB“Vk|2  (A.3) 

k=m+ 1 


Furthermore,  since  the  optimal  solution  fj  =  H  'b  allows  represen¬ 
tation  in  the  (dual)  form  fj  =  B§  (£  is  the  optimal  linear  combination 
of  the  representers),  the  upper  bound  of  the  terms  under  summation 
in  (A.3)  can  be  assessed  by 

\fjT B~^ <pk\  =  |fT0/<l  <  ffk(£TB  (A.4) 

Substituting  (A.4)  into  (A.3)  yields  the  following  upper  bound  on  the 
error  magnitude: 

em  <  f  TB~1£  J2  ak  ~  0(m~2)  (A.5) 

k>m 

This  estimate  remains  intact  if  we  assess  em  with  respect  to  the  norm 
induced  by  the  Hessian  matrix.  In  the  latter  case,  the  right-hand  side 
of  (A.5)  will  be  multiplied  by  a  scaling  factor  ||  H  ||/|  |B  '||  >  1. 


M.  Yaremchuket  al./ Ocean  Modelling97 (2016)  129-140 


139 


Dependence  of  the  distance  between  the  4dVar  solution  (Fig.  Alb) 
and  the  consecutive  approximations  to  the  a4dVar  solution  (Fig.  Aid) 
shown  in  the  inset  to  Fig.  A2,  confirms  the  above  estimate. 


Appendix  B.  H-orthogonalization  and  related  issues 


The  a4dVar  method  utilizes  the  technique  employed  by  Zupanski 
(2005)  in  the  Maximum  Likelihood  Ensemble  Filter,  which  is  based  on 
the  explicit  inversion  of  the  Hessian  matrix  in  the  subspace  spanned 
by  the  model  perturbations.  In  view  of  the  definition  (9),  B_1/2  can 
be  explicitly  represented  using  the  expression  for  the  square  root  of 
the  inverse  error  covariance: 

B~1/2  =  V  1  (I  -  yA)  (B.l ) 

which  allows  a  symmetric  Hessian  factorization 


H  =  HT/2H1/2 
where 


H1/2 


B-l/2  " 

Ho 

Hi  M1 


HjvM 


N 


(B.2) 


(B.3 ) 


is  the  Hessian  square  root. 

For  sufficiently  small  perturbations  8cm  =  epm,  perturbations  of 
the  auxiliary  vector 

SYm  =  H1/2Scm  (B.4) 


are  linear  in  Scm,  so  that  computation  of  the  dot  products  between 
the  vectors  8 Ym  provides  the  inner  product  in  the  control  space  as¬ 
sociated  with  the  Hessian  matrix 


SyJSY2  =  ScJhSc2  =  (<5ci ,  <5c2}p|,  (B.5) 

which  can  be  used  for  H-orthogonalization  of  the  search  subspaces 
of  the  a4dVar  algorithm. 

We  seek  the  optimal  correction  of  the  control  variable  c  in  the 
search  subspace  S  spanned  by  pm : 


for  each  m,  <5cJ  (b  -  He)  can  be  calculated  directly  from  the  variations 
of  the  cost  function  8Jm  =  J(c  +  Scm)  -J(c)  induced  by  <5cm: 

SJm  =  ^<5c^H<5cm  +  Scl (He  -  b) 

=  l3Yj«Ym-«4(b-Hc).  (B.9) 

Thus,  the  coefficients  for  the  optimal  correction  of  the  control 
variable  c  within  the  search  subspace  S  are  given  as  the  solution  to 
a  linear  system  posed  in  terms  of  the  quantities  8Jm  and  8Ym  com¬ 
puted  by  the  a4dVar  algorithm: 

n  , 

Y^SYlSY^^eUSYlSYrn-SJm).  (B.10) 

(=i  '  ' 

In  the  H-orthonormal  coordinate  system  5Y^i5Ym  =  e2,  and 
equations  (B.10)  are  simplified  to 

*  =  -t)'  (b"’ 

where  a/m  are  the  matrix  elements  of  the  linear  transformation  of  the 
original  basis  8cm  that  are  obtained  in  the  orthogonalization  process. 

For  a  differentiable  numerical  model  and  sufficiently  small  e,  the 
quadratic  term  in  the  right  hand  side  of  (B.10)  is  negligible.  In  the  re¬ 
ported  experiments  we  kept  it  in  place  since  the  value  of  e  was  close 
to  0.01  and  could  not  be  reduced  any  further  without  affecting  the 
rate  of  convergence.  The  relatively  large  limit  on  the  value  of  e  was 
caused  by  a  number  of  factors  deteriorating  the  linear  dependence 
between  the  magnitude  of  the  model  perturbations  and  e.  These  fac¬ 
tors  include  rounding  errors  (especially  for  temperature  and  salinity 
in  the  upper  layers),  non-differentiable  operators  in  the  model  code, 
particularly  at  the  open  boundary,  and  small-scale  instabilities  of  the 
flow  developing  at  the  mouth  of  the  Po  river  and  in  the  boundary 
currents,  especially  prominent  in  the  experiments  with  the  14-day 
assimilation  window. 

The  finite  value  of  e  also  affected  the  orthogonalization  process, 
resulting  in  non-zero  (of  the  order  of  3-10%)  off-diagonal  elements 
observed  in  the  Hessian  projection  after  orthogonalization  of  the  per¬ 
turbations.  For  that  reason,  the  optimal  coefficients  s(  were  computed 
through  the  direct  solution  of  Eq.  (B.10),  which  was  not  expensive  due 
to  the  low  dimension  of  the  system. 

References 


ms 

C  4-  C+  S/P/, 

/= 1 

where  the  coefficients  s(  satisfy  for  m  =  1.2 . ms. 

Pm  |h  (c  +  X)  S/P/)  -  bj  =0.  (B.6) 

This  constitutes  a  Ritz-Galerkin  projection  of  the  normal  system 
(4)  to  the  search  subspace,  S.  Rearranging,  we  obtain  the  linear  sys¬ 
tem  of  ms  equations  in  the  ms  unknowns  S} ,  s2, ,  sms: 

ms 

iSi  =  Pm(b-  He).  (B.7) 

1= 1 

Substituting  pm  =  Scm/s  into  (B.7),  multiplying  by  e2,  and  using 
(B.2)  and  (B.4)  yields 

m, 

^«Y^Y,S,  =  8«cJ(b-Hc).  (B.8) 

(=i 

The  right-hand  side  of  (B.8)  cannot  be  computed  directly  because 
evaluation  of  b  -  He  requires  the  adjoint  code  (Eq.  (5)).  Nonetheless, 


Ancell,  B.,  Hakim,  G.J.,  2007.  Comparing  adjoint-  and  ensemble-based  sensitivity  anal¬ 
ysis  with  applications  to  observation  targeting.  Mon.  Weather  Rev.  135, 4117-4134. 

Anderson,  J.L.,  Hoar,  T.,  Raeder,  K.,  Liu,  H.,  Collins,  N.,  Torn,  R.,  Arellano,  A.,  2009.  The 
data  assimilation  research  testbed:  a  community  facility .  Bull.  Am.  Meteorol.  Soc. 
90, 1283-1296. 

Barron,  C.N.,  Kara,  A.B.,  Hurlburt,  H.E.,  Rowley,  C.,  Smedstad,  L.F.,  2004.  Sea  surface 
height  predictions  from  the  global  Navy  Coastal  Ocean  Model  (NCOM)  during 
1998-2001.  J.  Atmos.  Oceanic  Technol.  21  (12),  1876-1894. 

Barron,  C.N.,  Kara,  A.B.,  Martin,  P.J.,  Rhodes,  R.C.,  Smedstad,  L.F.,  2006.  Formulation,  im¬ 
plementation  and  examination  of  vertical  coordinate  choices  in  the  global  Navy 
Coastal  Ocean  Model  (NCOM).  Ocean  Modell.  11,  347-375.  doi:10.1016/j.ocemod. 
2005.01.004. 

Bennett,  A.F.,  2002.  Inverse  Modeling  of  the  Ocean  and  Atmosphere.  Cambridge  Uni¬ 
versity  Press,  p.  234.  ISBN  0-521-81373-5. 

Buehner,  M.,  Houtekamer,  P.L.,  Charette,  C.,  Mitchell,  H.L.,  He,  B.,  2010.  Intercompari¬ 
son  of  variational  data  assimilation  and  the  ensemble  Kalman  filter  for  a  global 
deterministic  NWP.  Mon.  Wea.  Rev.  138, 1550-1566. 

Burrage,  D.M.,  Book,  J.W.,  Martin,  P.J.,  2009.  Eddies  and  filaments  of  the  Western  Adri¬ 
atic  cCurrent:  analysis  and  prediction.].  Mar.  Syst.  78,  S205-S226. 

Clayton,  A.M.,  Lorenc,  A.C.,  Barker,  D.M.,  2013.  Operational  implementation  of  a  hybrid 
ensemble/4D-Var  global  data  assimilation  system  at  the  met  office.  Q.  J.  R.  Meteo¬ 
rol.  Soc.  doi:10.1002/qj.2054. 

Cushman-Roisin,  B.,  Korotenko,  K.A.,  2007.  Mesoscale-resolving  simulations  of  summer 
and  winter  bora  events  in  the  Adriatic  Sea.  J.  Geophys.  Res.  112,  C11S91.  doi:  10. 

1029/2006JC003516. 

Desroziers,  G.,  Camino,  J.-T.,  Berre,  L.,  2014.  4DEnVar:  link  with  4D  state  formulation 
of  variational  assimilation  and  different  possible  implementations.  Q..  J.  R.  Meterol. 
Soc.  140,  2097-2110. 

Fairbairn,  D.,  Pring,  S.R.,  Lorenc,  A.C.,  Roulstone,  I.,  2014.  A  comparison  of  4dVar  with 
ensemble  data  assimilation  methods.  Q.  J.  R.  Meterol.  Soc.  140, 281-294. 


140 


M.  Yaremchuk  et  al.  /  Ocean  Modelling  97  (2016)  129-140 


Heimbach,  P.,  Hill,  C.,  Giering,  R.,  2005.  An  efficient  adjoint  of  the  parallel  MIT  GCM  gen¬ 
erated  via  automatic  differentiation.  Future  Gener.  Comput.  Syst.  21, 1356-1371. 

Holland,  W.R.,  Chow,  J.C.,  Bryan,  F.O.,  1998.  Application  of  a  third-order-upwind  scheme 
in  the  NCAR  Ocean  Model.  J.  Clim.  11, 1487-1493. 

Hoteit,  I.,  2008.  A  reduced-order  simulated  annealing  approach  for  four-dimensional 
variational  data  assimilation  in  meteorology  and  oceanography.  Int.  J.  Numer. 
Methods  Fluids  58, 1181-1199. 

Hoteit,  I.,  Hoar,  T.,  Gopalakrishnan,  G.,  Anderson,  J.,  Collins,  N.,  Cornuelle,  B.,  Kohl,  A., 
Heimbach,  P.,  2013.  A  MITgcm/DART  ocean  prediction  and  analysis  system  with 
application  to  the  Gulf  of  Mexico.  Dyn.  Atmos.  Oceans  63, 1-23. 

Hoteit,  I.,  Luo,  X.,  Pham,  D.T.,  2012.  Particle  Kalman  filtering:  a  nonlinear  Bayesian 
framework  for  ensemble  Kalman  filters.  Mon.  Weather  Rev.  140,  528-542. 

Ivatek-Sahdan,  S.,  Tudor,  M.,  2004.  Use  of  high-resolution  dynamical  adaptation  in  op¬ 
erational  suite  and  research  studies.  Meteorol.  Z.  13, 1-10. 

Kuhl,  D.D.,  Rosmond,  T.E.,  Bishop,  C.H.,  McLay,  J.,  Baker,  N.,  2013.  Comparison  of  hybrid 
ensemble/4DVar  and  4dVar  within  the  NAVDAS-AR  data  assimilation  framework. 
Mon.  Weather  Rev.  141,  2740-2758. 

Le  Dimet,  F.-X.,  Tlagrand,  O.,  1986.  Variational  algorithms  for  assimilation  and  analysis 
of  meteorological  observtaions:  theoretical  aspects.  Tellus  A  38A  (2),  97-110. 

Liu,  C.,  Xiao,  Q,,  Wang,  B.,  2008.  An  ensemble-based  four-dimensional  variational 
data  assimilation  scheme.  Part  I:  technical  formulation  and  preliminary  test.  Mon. 
Weather  Rev.  136,  3363-3373. 

Liu,  C.,  Xiao,  Q.,  Wang,  B.,  2009.  An  ensemble-based  four-dimensional  variational  data 
assimilation  scheme.  Part  II:  observing  system  simulation  experiments  with  ad¬ 
vanced  research  WRF  (ARW).  Mon. Weather  Rev.  137, 1687-1704. 

Martin,  P.J.,  2000.  A  Description  of  the  Navy  Coastal  Ocean  Model  Version  1.0.  NRL  Rep. 
NRL/FR/7322  00-9962, 42  pp.  Naval  Research  Laboratory,  Stennis  Space  Center,  MS. 

Martin,  P.J.,  Book,  J.W.,  Burrage,  D.M.,  Rowley,  C.D.,  Tudor,  M.,  2009.  Comparison  of 
model-simulated  and  observed  currents  in  the  central  adriatic  during  DART.  J.  Geo- 
phys.  Res.  114,  C01S05.  doi:10.1029/2008JC004842. 

Mellor,  G.L.,  1991.  An  equation  of  state  for  numerical  models  of  oceans  and  estuaries.  J. 
Atmos.  Oceanic  Technol.  8,  609-611. 

Mellor,  G.L.,  Yamada,  T.,  1974.  A  hierarchy  of  turbulence  closure  models  for  planetary 
boundary  layers.  J.  Atmos.  Sci.  31, 1791-1806. 

Menemenlis,  D.,  Wunsch,  G,  1997.  Linearization  of  an  oceanic  general  circulation 
model  for  data  assimilation  and  climate  studies.  J.  Atmos.  Oceanic  Technol.  14, 
1420-1443. 

Mirouze,  I.,  Weaver,  A.,  2010.  Representation  of  correlation  functions  in  variational  data 
assimilation  using  an  implicit  diffusion  operator.  Q,  J.  R.  Meteorol.  Soc.  136, 1421- 
1443. 

Moore,  A.M.,  Arango,  H.G.,  Broquet,  G.,  Powell,  B.S.,  Weaver,  A.T.,  Zavala-Garay,  J.,  2011. 
The  Regional  Ocean  Modeling  System  (ROMS)  4-dimensional  variational  data  as¬ 
similation  systems:  Part  I  System  overview  and  formulation.  Prog.  Oceanogr.  91 
(1),  34-49. 

Morey,  S.L.,  Martin,  P.J.,  OBrien,  J.J.,  Wallcraft,  A.A.,  Zavala-Hidalgo,  J.,  2003.  Export  path¬ 
ways  for  river  discharged  fresh  water  in  the  northern  Gulf  of  Mexico.  J.  Geophys. 
Res.  108  (CIO),  3303.  doi:10.1029/2002JC001674. 

Ngodock,  H.,  Carrier,  M.,  2014.  A  4dVar  system  for  the  navy  coastal  ocean  model.  Part 
I:  system  description  and  assimilation  of  synthetic  observations  in  Monterey  Bay. 
Mon.  Weather  Rev.  142  (6),  2085-2107. 

Oldenborgh,  G.J.,  Burgers,  G.,  Venzke,  S.,  Eckart,  G,  Giering,  R.,  1999.  Tracking  down  the 
ENSO  delayed  oscillator  with  an  adjoint  OGCM.  Mon.  Weather  Rev.  127, 1477-1495. 


Panteleev,  G.,  Yaremchuk,  M.,  Rogers,  E.,  2015.  Adjoint-free  variational  data  assimila¬ 
tion  into  a  regional  wave  model.  J.  Atmos.  Oceanic  Technol.  32, 1386-1399. 

Qui,  C.,  Shao,  A.,  Wei,  L.,  2007.  Fitting  model  fields  to  observations  by  using  singular 
value  decomposition:  an  ensemble-based  4dVar  approach.  J.  Geophys.  Res.  112, 
D11105.  doi:10.1029/2006JD007994. 

Reichel,  L.,  Ye,  Q.,  2005.  Breakdown-free  GMRES  for  singular  systems.  SIAM  J.  Matrix 
Anal.  Appl.  26  (4),  1001-1021. 

Robert,  C.,  Durbiano,  S.,  Blayo,  E.,  Verron,  J.,  Blum,  J.,  Dimet,  F.-X.  L.,  2005.  A  reduced- 
order  strategy  for  4dVar  data  assimilation.  J.  Mar.  Sys.  57  (1-2),  70-82. 

Rosmond,  T.,  Xu,  L.,  2006.  Development  of  the  NAVDAS-AR:  non-linear  formulation  and 
outer  loop  tests.  Tellus  58A,  45-58. 

Stammer,  D.,  Wunsch,  C.,  1996.  The  determination  of  the  large-scale  circulation  of  the 
Pacific  Ocean  from  satellite  altimetry  using  model  Green’s  functions.  J.  Geophys. 
Res.  101, 18409-18432. 

Tian,  X.,  Xie,  Z.,  2012.  Implementations  of  a  square-root  ensemble  analysis  and  a  hybrid 
localization  into  the  POD-based  ensemble  4DVar.  Tellus  A  64,  1-10.  doi:  10.3402/ 
tellusa.v64i0.18375. 

Torn,  R.D.,  Hakim,  G.J.,  2008.  Ensemble-based  sensitivity  analysis.  Mon.  Weather  Rev. 
136,  663-677. 

Trevisan,  A.,  DIsidoro,  M.,  Talagrand,  O.,  2010.  Four-dimensional  variational  assimila¬ 
tion  in  the  unstable  subspace  and  the  optimal  subspace  dimension.  Q,  J.  R.  Meteo- 
rol.  Soc.  136, 487-496. 

Weaver,  A.T.,  Courtier,  P.,  2001.  Correlation  modelling  on  a  sphere  using  a  generalized 
diffusion  equation.  Q.  J.  R.  Meteorol.  Soc.  127, 1815-1846. 

Weaver,  A.T.,  Deltel,  C.,  Machu,  E.,  Ricci,  S„  Daget,  N.,  2005.  A  multi-variate  balance 
operator  for  variational  data  assimilation.  Q.  J.  R.  Meteorol.  Soc.  131, 3605-3625. 

Weaver,  A.T.,  Mirouze,  I.,  2012.  On  the  diffusion  equation  and  its  application  to  isotropic 
and  anisotropic  correlation  modeling  in  variational  assimilation.  Q,  J.  R.  Meteorol. 
Soc.  138.  doi:10.1002/qj.l953. 

Xu,  L.,  Rosmond,  T.,  2004.  Formulation  of  the  NRL  Atmospheric  Variational  Data  Assim¬ 
ilation  System  -  Accelerated  Representer  (NAVDAS-AR).  Naval  research  Laboratory, 
p.  28.  NRL/MR/7532-36. 

Xu,  L.,  T.  Rosmond,  T.,  Daley,  R.,  2005.  Development  of  the  NAVDAS-AR:  formulation 
and  initial  tests  of  the  linear  problem.  Tellus  57A,  546-559. 

Yaremchuk,  M.,  Martin,  P.,  2014.  On  sensitivity  analysis  in  the  4dVar  framework.  Mon. 
Weather  Rev.  142,  774-787. 

Yaremchuk,  M.,  Nechaev,  D.,  2013.  Covariance  localization  with  the  diffusion-based 
correlation  models.  Mon.  Weather  Rev.  141,  848-860. 

Yaremchuk,  M.,  Nechaev,  D.,  Panteleev,  G.,  2009.  A  method  of  successive  corrections 
of  the  control  subspace  in  the  reduced-order  variational  data  assimilation.  Mon. 
Weather  Rev.  137, 2966-2978. 

Yaremchuk,  M.,  Sentchev,  A.,  2012.  Multi-scale  correlation  functions  associated  with 
polynomials  of  the  diffusion  operator.  Q,  J.  R.  Meterol.  Soc.  138, 1948-1953. 

Yaremchuk,  M.,  Smith,  S.,  2011.  On  the  correlation  functions  associated  with  polynomi¬ 
als  of  the  diffusion  operator.  Q.  J.  R.  Meterol.  Soc.  137, 1927-1932. 

Zhang,  F.,  Zhang,  M.,  Hansen,  J.A.,  2009.  Coupling  ensemble  Kalman  filter  with  four¬ 
dimensional  variational  data  assimilation.  Adv.  Atmos.  Sci.  26, 19. 

Zhang,  M.,  Zhang,  F.,  2012.  E4DVar:  Coupling  an  ensemble  Kalman  filter  with  four¬ 
dimensional  variational  data  assimilation  in  a  limited-area  weather  prediction 
model.  Mon.  Weather  Rev.  140,  587600. 

Zupanski,  M.,  2005.  Maximum  likelihood  ensemble  filter:  theoretical  aspects.  Mon. 
Weather  Rev.  133, 1710-1726. 


